This numerical exercise is based on Valeria Cristiani’s lecture. The goal of this numerical exercise is to familiarize ourselves with the different ways of identifying the components of a galaxy using galaxychop1:
JHistogram: Implementation of the dynamic decomposition model of galaxies described by Abadi et al.(2003).
JThreshold: Implementation of the dynamic decomposition model of galaxies used by Tissera et al.(2012), Vogelsberger et al.(2014), Marinacci et al.(2014), Park et al.(2019), etc.
KMeans: Implementation of Scikit-Learn K-means as a model for the dynamical decomposing of galaxies.
GaussianMixture: Implementation of the dynamic decomposition model of galaxies described by Obreja et al.(2018).
AutoGaussianMixture: Implementation of the dynamic decomposition model of galaxies described by Du et al.(2019).
import galaxychop as gchop # Installed and working
import numpy as np
import matplotlib.pylab as plt
import matplotlib.gridspec as gridspec
import matplotlib as mpl
import pandas as pd
%matplotlib inline
plt.style.use('dark_background')
plt.style.use('seaborn-bright')
COLOR = 'lightgray'
mpl.rcParams['text.color'] = COLOR # labels, titles, etc colors
mpl.rcParams['axes.labelcolor'] = COLOR # Numbers in axis colors
mpl.rcParams['xtick.color'] = COLOR # ticks colors
mpl.rcParams['ytick.color'] = COLOR # ticks colors
mpl.rc( 'axes',edgecolor = COLOR) # axis Colors
mpl.rc('font', size=15) # Font Size
mpl.rcParams['figure.figsize'] = [10, 10] # Size of figures
gal = gchop.read_hdf5("gal394242.h5")
gal2 = gchop.read_hdf5("gal443801.h5")
print(f'SubhaloID 394242: {gal}\nSubhaloID 443801: {gal2}\n')
SubhaloID 394242: <Galaxy stars=37393, dark_matter=155101, gas=80153, potential=True> SubhaloID 443801: <Galaxy stars=47789, dark_matter=282694, gas=4459, potential=True>
if not gchop.preproc.is_centered_and_aligned(gal):
print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal))
#gal = gchop.preproc.center_and_align(gal)
gal = gchop.preproc.center_and_align(gal, r_cut=150)
print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal))
Centered & Aligned: False Centered & Aligned: True
if not gchop.preproc.is_centered_and_aligned(gal2):
print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal2))
#gal2 = gchop.preproc.center_and_align(gal2)
gal2 = gchop.preproc.center_and_align(gal2, r_cut=260)
print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal2))
Centered & Aligned: False Centered & Aligned: True
In the jupyter lab I have run:
########################################
# Packages
import numpy as np
import matplotlib.pylab as plt
import illustris_python as il
import matplotlib.gridspec as gridspec
########################################
# Getting the data
basePath = '/home/tnguser/sims.TNG/TNG100-1/output/'
snapNum = 99 # z=0
fields = ['SubhaloFlag', 'SubhaloHalfmassRadType',
'SubhaloMassType', 'SubhaloPos', 'SubhaloVel',
'SubhaloSFRinRad', 'SubhaloSFR','SubhaloStarMetallicity','SubhaloGrNr']
Subhalo = il.groupcat.loadSubhalos(basePath=basePath,snapNum=snapNum,fields=fields)
FIELDS=['GroupFirstSub','GroupMassType','Group_M_Crit200','Group_R_Crit200']
HHH=il.groupcat.loadHalos(basePath,99,fields=FIELDS)
GFSub = HHH['GroupFirstSub']
Mtype = HHH['GroupMassType']
M200 = HHH['Group_M_Crit200']
R200 = HHH['Group_R_Crit200']
########################################
# Saving the data
gdms=np.array([0,1,4]) # to keep just columns of gas,dm and stars
Data={}
for i in [394242, 443801]:
a=Subhalo['SubhaloGrNr'][i]
Data[str(i)]= dict( id = i, # The subhalo "ID" is the same as the "index"
cpos = Subhalo['SubhaloPos'][i],
cvel = Subhalo['SubhaloVel'][i],
rh_type = Subhalo['SubhaloHalfmassRadType'][i][gdms],
M_type = Subhalo['SubhaloMassType'][i][gdms],
sfr_2rh = Subhalo['SubhaloSFRinRad'][i], # in 2 * R_half
sfr_tot = Subhalo['SubhaloSFR'][i], # in 2 * R_half
mh = Subhalo['SubhaloStarMetallicity'][i],
Central = (i in GFSub),
M200_parent = M200[a],
R200_parent = R200[a])
np.save('Gal_Data.npy', Data)
Now from this file we can print the info
Data=np.load('Gal_Data.npy',allow_pickle=True).item()
pd.DataFrame(Data)
| 394242 | 443801 | |
|---|---|---|
| id | 394242 | 443801 |
| cpos | [63794.09, 12026.507, 5045.7017] | [15467.963, 63014.26, 73353.84] |
| cvel | [194.05467, 263.4645, 114.92622] | [157.85007, -136.74078, -32.544792] |
| rh_type | [60.96873, 65.7992, 5.139353] | [109.72671, 74.55964, 1.619801] |
| M_type | [8.251568, 78.41508, 2.5353656] | [0.48876184, 142.92282, 3.0524666] |
| sfr_2rh | 0.355641 | 0.0 |
| sfr_tot | 0.93635 | 0.0 |
| mh | 0.01982 | 0.024343 |
| Central | False | True |
| M200_parent | 190.210754 | 139.27948 |
| R200_parent | 201.506332 | 181.628571 |
x, y, z = gal.stars.x.value , gal.stars.y.value , gal.stars.z.value
r= (x**2+y**2+z**2)**0.5
# print(sum(r>70))
gal.stars.to_dataframe()
| ptype | ptypev | m | x | y | z | vx | vy | vz | softening | potential | kinetic_energy | total_energy | Jx | Jy | Jz | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | stars | 0 | 5.224283e+05 | -0.046012 | -0.062069 | -0.262818 | 6.362919 | 22.511274 | -5.774254 | 0.0 | -195699.620206 | 290.293111 | -195409.327095 | 6.274766 | -1.937972 | -0.640845 |
| 1 | stars | 0 | 9.745897e+05 | -0.155791 | -0.125822 | -0.183540 | 21.987698 | 8.126047 | 0.272632 | 0.0 | -196176.962277 | 274.782915 | -195902.179362 | 1.457155 | -3.993156 | 1.500567 |
| 2 | stars | 0 | 6.935776e+05 | -0.250341 | -0.162899 | -0.206341 | -9.982170 | 7.308657 | 15.308136 | 0.0 | -195152.120168 | 193.699612 | -194958.420557 | -0.985599 | 5.891984 | -3.455738 |
| 3 | stars | 0 | 1.070959e+06 | -0.223984 | -0.253224 | -0.378401 | -9.688729 | -5.382013 | -1.178408 | 0.0 | -194695.767625 | 62.113089 | -194633.654536 | -1.738158 | 3.402281 | -1.247935 |
| 4 | stars | 0 | 6.013803e+05 | -0.068238 | -0.198522 | -0.285535 | 17.203435 | 7.988667 | -16.385347 | 0.0 | -195703.942688 | 314.128285 | -195389.814402 | 5.533893 | -6.030286 | 2.870127 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 37388 | stars | 0 | 9.867470e+05 | 134.949071 | -161.114980 | -11.970927 | -59.939303 | 57.887495 | 120.366415 | 0.0 | -26553.893402 | 10715.877926 | -15838.015476 | -18699.865471 | -15525.806786 | -1845.255888 |
| 37389 | stars | 0 | 8.290853e+05 | 15.012682 | 15.493519 | -4.697090 | -13.645752 | 262.439571 | 314.563208 | 0.0 | -101280.155757 | 84005.373302 | -17274.782455 | 6106.393180 | -4658.342086 | 4151.342534 |
| 37390 | stars | 0 | 1.004103e+06 | 17.379644 | 56.147464 | -2.132304 | -41.672564 | 304.593269 | 11.333777 | 0.0 | -64131.846872 | 47321.058282 | -16810.788589 | 1285.848336 | -108.118416 | 7633.531240 |
| 37391 | stars | 0 | 1.526682e+06 | 269.301518 | -82.400695 | 23.001600 | -88.927482 | -50.437606 | -26.101538 | 0.0 | -19698.409864 | 5566.669709 | -14131.740154 | 3310.930548 | 4983.709501 | -20910.610228 |
| 37392 | stars | 0 | 1.356703e+06 | 13.433154 | -8.906366 | 3.706972 | 172.602906 | 380.189578 | -150.006406 | 0.0 | -111482.516414 | 98418.900084 | -13063.616330 | -73.340046 | 2654.893305 | 6644.409980 |
37393 rows × 16 columns
PLOT=True
LIM=6
if PLOT: # High Resolution
x, y, z = gal.stars.x.value , gal.stars.y.value , gal.stars.z.value
mask = (((x**2+y**2+z**2)**0.5)<LIM)
x, y, z = x[mask], y[mask] , z[mask]
H, edges=np.histogramdd((x,y,z),bins=50) #bins=(50,45,55)
x= (edges[0][1:] + edges[0][:-1])/2
y= (edges[1][1:] + edges[1][:-1])/2
z= (edges[2][1:] + edges[2][:-1])/2
x,y,z = np.meshgrid(x,y,z)
x,y,z=x.ravel(),y.ravel(),z.ravel()
#------------------------------------------
x2, y2, z2 = gal2.stars.x.value , gal2.stars.y.value , gal2.stars.z.value
mask = (((x2**2+y2**2+z2**2)**0.5)<LIM)
x2, y2, z2 = x2[mask], y2[mask] , z2[mask]
H2, edges=np.histogramdd((x2,y2,z2),bins=50) #bins=(50,45,55)
x2= (edges[0][1:] + edges[0][:-1])/2
y2= (edges[1][1:] + edges[1][:-1])/2
z2= (edges[2][1:] + edges[2][:-1])/2
x2,y2,z2 = np.meshgrid(x2,y2,z2)
x2,y2,z2=x2.ravel(),y2.ravel(),z2.ravel()
if PLOT: # Scatter density version (plotly)
import plotly.graph_objects as go
from plotly.subplots import make_subplots
#------------------
# Figure
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'volume'},{'type': 'volume'}]],
print_grid=False,
subplot_titles=( ['SubHalo: 394242','SubHalo: 443801']),
horizontal_spacing = 0.075
)
#------------------
# Scatter plot High resolution
trace1 = go.Volume(name = '394242',
x = y,
y = x,
z = z,
value = H.ravel(),
isomin = 1,
isomax = H.max(),
opacity = 0.15, # needs to be small to see through all surfaces
opacityscale = 'uniform',#""max",
surface_count= 25, # needs to be a large number for good volume rendering
colorscale = 'matter',
showscale = False,
# colorbar = dict(thickness = 20,
# outlinewidth = 0.7,
# title = "Count particles per bin",
# titleside = "right",
# tickmode = "array",
# tickvals = [], #[0.001,0.999]
# #labelalias = {1: str(H.max()), 0:str(H.min()) },
# ticks = "outside",
# x = -0.2,) # position of colorbar
)
#------------------
# Scatter plot High resolution
trace2 = go.Volume(name = '443801',
x = y2,
y = x2,
z = z2,
value = H2.ravel(),
isomin = 1,
isomax = H2.max(),
opacity = 0.15, # needs to be small to see through all surfaces
opacityscale = 'uniform',#"max",
surface_count= 25, # needs to be a large number for good volume rendering
colorscale = 'matter',
colorbar = dict(thickness = 20,
outlinewidth = 0.7,
title = "Count Stars",
titleside = "right",
tickmode = "array",
tickvals = [], #[0.001,0.999]
#labelalias = {1: str(H.max()), 0:str(H.min()) },
ticks = "outside",
x = 1.2,) # position of colorbar
)
fig.add_trace(trace1, # trace= : figure
row=1, # row= : Position y
col=1) # col= : Position x
fig.add_trace(trace2, # trace= : figure
row=1, # row= : Position y
col=2) # col= : Position x
#------------------
# Background color (page and axes)
# https://www.w3schools.com/colors/colors_rgb.asp
fig['layout']['paper_bgcolor']="rgb(17,17,17)"
set_bgcolor= dict(showbackground = True,
backgroundcolor = "rgb(101, 101, 101)", #"rgb(212, 219, 249)"
gridcolor = "rgb(10, 10, 10)",
zeroline = False)
fig.update_scenes(xaxis=set_bgcolor,
yaxis=set_bgcolor,
zaxis=set_bgcolor)
#------------------
# Font color and size
fig.update_layout(font_color = 'white', font_size = 16)
fig.update_layout(height=550, width=1000,title_text="Stars in Halos")
fig.show()
print()